!     -------------
!     PROGRAM JANUS   GF-9812
!     -------------
!
! --- Object -----------------------------------------------------------
!
!     This program computes the lunar coordinates with the tables and 
!     algorithms given in :
!     Lunar Tables and Programs from 4000 B.C. to A.D. 8000
!     M. Chapront-Touze et J. Chapront (Willmann Bell, 1991).
!
!     The lunar tables are stored in the file : janus.dat.
!     The results are stored in the file : janus.txt.
!
!V98  REMARK:
!V98  THE VERSION 98-12 IS DIFFERENT FROM THE ORIGINAl VERSION :
!V98  THERE IS NO MORE DISTINCTION BETWEEN 3 LEVELS OF PRECISION.
!V98  THE HIGHEST PRECISION IS THE ONLY ONE USED WHEN READING TABLES.
!V98  THESE CHANGES ARE NOTED BY 'V98'. 
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
      character rep
      character*30 datcal
!
! --- Initialization ---------------------------------------------------
!
      data d2000/2451545.d0/
      data sjul/36525.d0/
!
      nultab=4
      nulprn=10
!
      ipr=1
      icd=1
      ir=1
      idt=1
      idf=2
      d=2451545.0d0
      datcal='01 Janv 2000   12h00m00s'
      ndat=10
      pas=1.0d0
!
! --- Results File -----------------------------------------------------
!
      open (nulprn,file='janus.txt')
      write (nulprn,1000)
!
! --- Reading the Lunar Tables -----------------------------------------
!
      call LECTAB (1,1,nultab)
!
! --- Reading the parameters of computation ----------------------------
!
100   continue
!
      call MENU (ipr,icd,ir,idt,idf,d,datcal,ndat,pas)
!V98  call LECTAB (2,ipr,nultab)
!
! --- Computation : Start ----------------------------------------------
!
      t0=d
      dt=0.d0
      pasj=pas/24.d0
!
      select case (icd)
!
! --- Compute the lunar coordinates ------------------------------------
!
      case (1,2)
!
      do i=1,ndat
!
         if (mod(i,10).eq.1) call TITRE (1,nulprn,ipr,icd,ir,idt,idf)
         if (i.eq.1) call TITRE (2,nulprn,ipr,icd,ir,idt,idf)
!
         dj=t0+(i-1)*pasj
         t=(dj-d2000)/sjul
!
         if (idt.eq.2) then
            call UT (dj,dt,ierr,limut,tlimut)
            t=t+dt/315576.d4
         endif
!
         call LUNPOL (t,c1,c2,c3,ir,icd,ipr)
!
         select case (ir)
         case (1)
            call IMPC (dj,nulprn,c1,c2,c3,ipr,dt,idt,idf)
         case (2)
            call IMPQ (dj,nulprn,c1,c2,c3,ipr,dt,idt,idf)
         case (3)
            call IMPC (dj,nulprn,c1,c2,c3,ipr,dt,idt,idf)
         case (4)
            call IMPQ (dj,nulprn,c1,c2,c3,ipr,dt,idt,idf)
         endselect
!
         if (mod(i,10).eq.0) then
            write (*,'(/10x,a,$)') ' Press Enter '
            read (*,*)
         endif
!
      enddo
!
! --- Compute the osculating elements of the lunar orbit ---------------
!
      case (3)
!
      do i=1,ndat
!
         if (mod(i,10).eq.1) call TITRE (1,nulprn,ipr,icd,ir,idt,idf)
         if (i.eq.1) call TITRE (2,nulprn,ipr,icd,ir,idt,idf)
!
         dj=t0+(i-1)*pasj
         t=(dj-d2000)/sjul
!
         call LUNOSC (t,a,e,gam,pom,gom,alam)
         call IMPS (dj,nulprn,a,e,gam,pom,gom,alam,ipr,idf)
!
         if (mod(i,10).eq.0) then
            write (*,'(/10x,a,$)') ' Press Enter '
            read (*,*)
         endif
!
      enddo
!
! --- Computation : End ------------------------------------------------
!
      endselect
!
! --- Another computation ? --------------------------------------------
!
      write (*,1001)
      read (*,1002) rep
!
      if (rep.eq.'n'.or.rep.eq.'N') then
         write (nulprn,'(//)')
         stop
      endif
!
      goto 100
!
! --- Formats ----------------------------------------------------------
!
1000  format (2x,'JANUS GF9812',15x,'MOON from -4000 to 8000',16x,
     .        'M&J CHAPRONT'/2x,78('-'))
1001  format (/6x,'Done, Another one (y/n) ? ',$)
1002  format (a)
!
      end
!
!
!
      BLOCK DATA
! ----------------------------------------------------------------------
!
!     Ref : GF-9812
!
!     Secular terms of the mean longitude of the Moon.
!
! ----------------------------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
      common/long/al(5)
!
      data al/
     .218.31665d0,481267.88134d0,-13.268d-4,1.856d-6,-1.534d-8/
!
      end
!
!
!
      subroutine MENU (ipr,icd,ir,idt,idf,d,datcal,ndat,pas)
!-----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Object -----------------------------------------------------------
!
!     Display the Menu.
!
! --- Input ------------------------------------------------------------
!
!     ipr      Index of Precision (integer).
!     icd      Index of Coordinates (integer).
!     ir       Index of Frame (integer).
!     idt      Index of Time (integer).
!     idf      Index of Date (integer).
!     d        Initial Julian date (real*8).
!     datcal   Inititial Calendar date (character*30).
!     ndat     Number of dates (integer).
!     pas      Time interval in hours (real*8).
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
      character rep
      character*15 z
      character*30 apr(3),acd(3),arp(4),adt(2),datcal
!
! --- Initialization ---------------------------------------------------
!
      data apr/
     .'High precision',
     .'Middle precision',
     .'Low precision'/
      data acd/
     .'Geometric Coordinates',
     .'Apparent Coordinates',
     .'Osculating Elements'/
      data arp/
     .'Mean Ecliptic Frame',
     .'Mean Equatorial Frame',
     .'True Ecliptic Frame',
     .'True Equatorial Frame'/
      data adt/
     .'Dynamical Time (ET)',
     .'Universal Time (UT)'/
!
! --- Display the Menu -------------------------------------------------
!
100   continue
!
      call EFFECR
      write (*,1000)
!
!V98  write (*,1001) apr(ipr)
      write (*,1002) acd(icd)
      write (*,1003) arp(ir)
      write (*,1004) adt(idt)
!
      if (idf.eq.1) write (*,1005) d
      if (idf.eq.2) write (*,1006) datcal
!
      write (z,'(i15)') ndat
      i=1
      do while (z(i:i).eq.' ')
         i=i+1
      enddo
      write (*,1007) z(i:15)
!
      write (z,'(f15.4)') pas
      i=1
      do while (z(i:i).eq.' ')
         i=i+1
      enddo
      if (z(i:i).eq.'.') then
         i=i-1
         z(i:i)='0'
      endif
      write (*,1008) z(i:15)
!
      ierr=0
      write (*,1009)
      read (*,'(a)') rep
!
      select case (rep)
!V98  case ('P','p')
!V98     call LECPRM (1,ipr,icd,ir,idt,idf,d,datcal,ndat,pas)
      case ('C','c')
         call LECPRM (2,ipr,icd,ir,idt,idf,d,datcal,ndat,pas)
      case ('F','f')
         if (icd.eq.3) then
            ir=1
            write (*,2001)
            read (*,*)
         else
            call LECPRM (3,ipr,icd,ir,idt,idf,d,datcal,ndat,pas)
         endif
      case ('T','t')
         if (icd.eq.3) then
            idt=1
            write (*,2002)
            read (*,*)
         else
            call LECPRM (4,ipr,icd,ir,idt,idf,d,datcal,ndat,pas)
         endif
      case ('D','d')
         call LECPRM (5,ipr,icd,ir,idt,idf,d,datcal,ndat,pas)
      case ('N','n')
         call LECPRM (6,ipr,icd,ir,idt,idf,d,datcal,ndat,pas)
      case ('I','i')
         call LECPRM (7,ipr,icd,ir,idt,idf,d,datcal,ndat,pas)
      case ('X','x')
         return
      case ('Q','q')
         stop
      case default
         ierr=1
      endselect
!
      if (ierr.eq.1) then
         write (*,2003)
         read (*,*)
      endif
!
      goto 100
!
! --- Formats ----------------------------------------------------------
!
1000  format (2x,'JANUS GF9812',15x,'MOON from -4000 to 8000',16x,
     .        'M&J CHAPRONT'/2x,78('-')//
     .        8x,'LIST OF COMMANDS',t40,'PARAMETERS OF COMPUTATION'/)
!V98  1001 format (8x,'P : Precision       ',t40,a)
1002  format (8x,'C : Coordinates          ',t40,a)
1003  format (8x,'F : Frame                ',t40,a)
1004  format (8x,'T : Time                 ',t40,a)
1005  format (8x,'D : Initial Date         ',t40,'JD',f13.5)
1006  format (8x,'D : Initial Date         ',t40,a)
1007  format (8x,'N : Number of dates      ',t40,a)
1008  format (8x,'I : Time Interval (hour) ',t40,a)
1009  format (8x,'X : Run                  '/
     .        8x,'Q : Quit                 '//
     .        8x,'Enter your choice ? ',$)
!
2001  format (/8x,'For the osculating elements : ',
     .        'Mean Ecliptic of the date'/
     .        /8x,'Press Enter ',$)
2002  format (/8x,'For the osculating elements : ',
     .        'Dynamical Time (ET)'/
     .        /8x,'Press Enter ',$)
2003  format (/8x,'Error : Press Enter ',$)
!
      end
!
!
!
      subroutine LECTAB (ilec,ipr,nul)
!-----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Object -----------------------------------------------------------
!
!     Reading the lunar tables : Longitude, Latitude and Distance.
!
! --- Input ------------------------------------------------------------
!
!     ilec     Index of reading (integer).
!              ilec = 1 : initial reading of the lunar tables.
!              ilec = 2 : reading according to the level of precision.
!
!     ipr      Index of precision (integer).
!              ipr = 1 High precision.
!              ipr = 2 Middle precision.
!              ipr = 3 Low precision.
!
!     nul      Logical unit number of the file 'lunar tables' (integer).
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
      character*80 zlec(1324)
!
      dimension nc(3,3),np(3,3,3)
      dimension a(5),b(2)
!
      common/cv/deg,hr
      common/lu/pbc(218,3),apbc(218,5,3),per(423,3),aper(423,2,3),
     .          kc(3),kp(3,3)
!
! --- Initialization ---------------------------------------------------
!
      data pi/3.141592653589793d0/
!
      data nc/218,188,154,59,45,40,29,14,18/
!
      data np/
     .244,64,114,154,64,68,25,12,19,
     .3,2,0,16,6,9,5,0,4,1,0,0,6,0,3,1,0,0/
!
      deg=pi/180.d0
      hr=deg*15
      il=0
!
! --- Initial reading of the lunar tables ------------------------------
!
      if (ilec.eq.1) then
         call EFFECR
         write (*,1000)
         write (*,1001)
         open (nul,file='janus.dat',status='old')
         do i=1,1324
            read (nul,2000) zlec(i)
         enddo
         close (nul)
         write (*,1002)
         read (*,*)
      endif
!
! --- Reading according to the level of precision -----------------------
!
      do i=1,3
!
         kc(i)=nc(i,ipr)
         jmax=kc(i)
!
         do j=1,jmax
!
            il=il+1
            if (i.ne.3) then
               read (zlec(il),2001) c,(a(j1),j1=1,5)
            else
               read (zlec(il),2002) c,(a(j1),j1=1,5)
               a(1)=a(1)+90.d0
            endif
!
            pbc(j,i)=c
!
            do j1=1,2
               apbc(j,j1,i)=a(j1)*deg
            enddo
!
            apbc(j,3,i)=a(3)*deg*1.d-4
            apbc(j,4,i)=a(4)*deg*1.d-6
            apbc(j,5,i)=a(5)*deg*1.d-8
!
         enddo
!
         jsup=nc(i,1)-jmax
!
         do j=1,jsup
            il=il+1
            read (zlec(il),2003)
         enddo
!
         jmax=0
!
         do k=1,3
!
            kp(i,k)=np(i,k,ipr)
            jmin=jmax+1
            jmax=jmax+kp(i,k)
!
            do j=jmin,jmax
!
               il=il+1
               if (i.ne.3) then
                  if (k.eq.1) read (zlec(il),3001) c,(b(j1),j1=1,2)
                  if (k.eq.2) read (zlec(il),4001) c,(b(j1),j1=1,2)
                  if (k.eq.3) then
                     read (zlec(il),5001) c,(b(j1),j1=1,2)
                     c=c*1.d-4
                  endif
                  c=c*1.d-3
               else
                  if (k.eq.1) read (zlec(il),3002) c,(b(j1),j1=1,2)
                  if (k.eq.2) read (zlec(il),4002) c,(b(j1),j1=1,2)
                  if (k.eq.3) then
                     read (zlec(il),5002) c,(b(j1),j1=1,2)
                     c=c*1.d-4
                  endif
                  b(1)=b(1)+90
               endif
!
               per(j,i)=c
!
               do j1=1,2
                  aper(j,j1,i)=b(j1)*deg
               enddo
!
            enddo
!
            jsup=np(i,k,1)-kp(i,k)
!
            do j=1,jsup
               il=il+1
               read (zlec(il),2003)
            enddo
!
         enddo
!
      enddo
!
      return
!
! --- Formats ----------------------------------------------------------
!
1000  format (///
     .        2x,t20,'          PROGRAM JANUS GF9812        '//
     .        2x,t20,'               REFERENCE:             '//
     .        2x,t20,'       LUNAR TABLES AND PROGRAMS      '/
     .        2x,t20,'      FROM 4000 B.C. TO A.D. 8000     '//
     .        2x,t20,'   M. CHAPRONT-TOUZE AND J. CHAPRONT  '//)
1001  format (2x,t20,'      Reading the lunar tables ...    '//)
1002  format (2x,t20,'              Press Enter             ',$)
!
2000  format (a)
2001  format(5x,f11.8,2x,f11.7,2x,f16.7,2(2x,f8.3),2x,f7.3)
2002  format(5x,f11.4,2x,f11.7,2x,f16.7,2(2x,f8.3),2x,f7.3)
2003  format(1x)
!
3001  format(5x,f8.5,2x,f8.4,2x,f13.4)
3002  format(5x,f8.4,2x,f8.4,2x,f13.4)
!
4001  format(5x,f8.5,2x,f7.3,2x,f12.3)
4002  format(5x,f8.4,2x,f7.3,2x,f12.3)
!
5001  format(5x,f5.2,2x,f5.1,2x,f10.1)
5002  format(5x,f4.1,2x,f5.1,2x,f10.1)
!
      end
!
!
!
      subroutine LECPRM (iprm,ipr,icd,ir,idt,idf,d,datcal,ndat,pas)
!-----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Object -----------------------------------------------------------
!
!     Reading the parameters of the computation.
!
! --- Input ------------------------------------------------------------
!
!     iprm     Index of the parameters (integer).
!              1 : ipr    Precision.
!              2 : icd    Coordinates.
!              3 : ir     Frame.
!              4 : idt    Time.
!              5 : idf    Julian date / Calendar date.
!              5 : d      Initial Julian date.
!              5 : datcal Initial Calendar date.
!              6 : ndat   Number of dates.
!              7 : pas    Time interval in hours.
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
      character*16 adat,apas
      character*30 datcal,datx
!
! --- Initialization ---------------------------------------------------
!
      djmin=260057.0d0
      djmax=4643366.d0
      datmin=-40000101.0d0
      datmax=80010101.0d0
      call UT (d,dt,ierr,limut,tlimut)
!
! --- Select the parameters : Start ------------------------------------
!
      select case (iprm)
!
! --- Precision --------------------------------------------------------
!V98
!V98  case (1)
!V98
!V98  write (*,1001)
!V98  read (*,'(i10)',iostat=nerr) ix
!V98  if (nerr.ne.0) goto 100
!V98  if (ix.lt.1.or.ix.gt.3) goto 100
!V98  ipr=ix
!V98!
! --- Coordinates ------------------------------------------------------
!
      case (2)
!
      write (*,1002)
      read (*,'(i10)',iostat=nerr) ix
      if (nerr.ne.0) goto 100
      if (ix.lt.1.or.ix.gt.3) goto 100
      icd=ix
      if (icd.eq.3) then
         ir=1
         idt=1
      endif
!
! --- Frame ------------------------------------------------------------
!
      case (3)
!
      write (*,1003)
      read (*,'(i10)',iostat=nerr) ix
      if (nerr.ne.0) goto 100
      if (ix.lt.1.or.ix.gt.4) goto 100
      ir=ix
!
! --- Time -------------------------------------------------------------
!
      case (4)
!
      write (*,1004)
      read (*,'(i10)',iostat=nerr) ix
      if (nerr.ne.0) goto 100
      if (ix.lt.1.or.ix.gt.2) goto 100
!
      df=d+pas*(ndat-1)/24.d0
      if (ix.eq.2.and.(df-tlimut).gt.0.d0) goto 400
      idt=ix
!
! --- Initial date -----------------------------------------------------
!
      case (5)
!
      write (*,1005)
      read (*,'(i10)',iostat=nerr) ix
      if (nerr.ne.0) goto 100
      if (ix.lt.1.or.ix.gt.2) goto 100
      idf=ix
!
      if (idf.eq.1) then
         write (*,1006)
         read (*,'(f20.8)',iostat=nerr) x
         if (nerr.ne.0) goto 100
         if (x.eq.0.d0) return
         if (x.lt.djmin.or.x.gt.djmax) goto 100
      else
         write (*,1007)
         read (*,'(a)',iostat=nerr) adat
         if (nerr.ne.0) goto 100
         if (adat.eq.' ') return
         n=index(adat,'.')
         if (n.eq.0) then
            read (adat,'(i16)',iostat=nerr) ix
            if (nerr.ne.0) goto 100
            xx=ix
         else
            read (adat,'(f16.6)',iostat=nerr) xx
            if (nerr.ne.0) goto 100
         endif
         if (xx.lt.datmin.or.xx.gt.datmax) goto 100
         call DCDJ (xx,x,ierr,datx)
         if (ierr.eq.1) goto 100
         if (ierr.eq.2) goto 200
      endif
!
      df=x+pas*(ndat-1)/24.d0
      if ((df-djmax).gt.0.d0) goto 300
      if (idt.eq.2.and.(df-tlimut).gt.0.d0) goto 400
      d=x
      call DJDATE (d,3,ij,im,ia,nh,nm,ns)
      date=abs(ia)*1.d4+im*1.d2+ij+nh*1.d-2+nm*1.d-4+ns*1.d-6
      if (ia.lt.0) date=-date
      call DCDJ (date,x,ierr,datcal)
!
! --- Number of dates --------------------------------------------------
!
      case (6)
!
      write (*,1008)
      read (*,'(i10)',iostat=nerr) ix
      if (nerr.ne.0) goto 100
      if (ix.lt.1.or.ix.gt.1000) goto 100
!
      df=d+pas*(ix-1)/24.d0
      if ((df-djmax).gt.0.d0) goto 300
      if (idt.eq.2.and.(df-tlimut).gt.0.d0) goto 400
      ndat=ix
!
! --- Time interval ----------------------------------------------------
!
      case (7)
!
      write (*,1009)
      read (*,'(a)',iostat=nerr) apas
      if (ierr.ne.0) goto 100
      if (apas.eq.' ') goto 100
      n=index(apas,'.')
      if (n.eq.0) then
         read (apas,'(i16)',iostat=nerr) ix
         if (nerr.ne.0) goto 100
         x=ix
      else
         read (apas,'(f16.8)',iostat=nerr) x
         if (nerr.ne.0) goto 100
      endif
      if (x.le.0.d0.or.x.gt.1000.d0) goto 100
!
      df=d+x*(ndat-1)/24.d0
      if ((df-djmax).gt.0.d0) goto 300
      if (idt.eq.2.and.(df-tlimut).gt.0.d0) goto 400
      pas=x
!
! --- Select the parameters : End --------------------------------------
!
      endselect
      return
!
! --- Errors Messages --------------------------------------------------
!
100   continue
!
      write (*,2001)
      read (*,*)
      return
!
200   continue
!
      write (*,2002)
      read (*,*)
      return
!
300   continue
!
      write (*,2003)
      read (*,*)
      return
!
400   continue
!
      write (*,2004) limut
      read (*,*)
      return
!
! --- Formats ----------------------------------------------------------
!
!V98  1001  format (/8x,'1: High precision'
!V98       .        /8x,'2: Middle precision'
!V98       .        /8x,'3: Low precision',
!V98       .        t40,'Your choice ? ',$)
1002  format (/8x,'1: Geometric Polar Coordinates'
     .        /8x,'2: Apparent  Polar Coordinates'
     .        /8x,'3: Osculating Elements',
     .        t50,'Your choice ? ',$)
1003  format (/8x,'1: Mean Ecliptic and Equinox of the date'
     .        /8x,'2: Mean Equator  and Equinox of the date'
     .        /8x,'3: Mean Ecliptic of the date and True Equinox'
     .        /8x,'4: True Equator  and True Equinox',
     .        t60,'Your choice ? ',$)
1004  format (/8x,'1: Dynamical Time (ET)'
     .        /8x,'2: Universal Time (UT)',
     .        t40,'Your choice ? ',$)
1005  format (/8x,'1: Julian   Date'
     .        /8x,'2: Calendar Date ',
     .        t40,'Your choice ? ',$)
1006  format (/8x,'Enter a Julian date (decimal) : ',$)
1007  format (/8x,'Enter a Calendar date ',
     .        '(+/-yyyymmdd.hhmmss) : ',$)
1008  format (/8x,'Enter the number of dates (integer) : ',$)
1009  format (/8x,'Enter the time interval en hours (decimal) : ',$)
!
2001  format (/8x,'Error : Press Enter ',$)
2002  format (/8x,'Error : No valid date in Gregorian calendar',
     .        ' : Press Enter ',$)
2003  format (/8x,'Error : Date after 8000',
     .        ' : Press Enter ',$)
2004  format (/8x,'Erreur : Date after UT limit',i4,
     .        ' : Press Enter ',$)
!
      end
!
!
!
      subroutine LUNPOL (t,c1,c2,c3,ir,icd,ipr)
! ----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Object -----------------------------------------------------------
!
!     Polar Coordinates of the Moon.                
!
! --- Input ------------------------------------------------------------
!
!     t        Time in Julian century since J2000 (real*8).
!
!     ir       Index of Frame (integer).
!              ir = 1 : Mean Ecliptic and Equinox of the date.
!              ir = 2 : Mean Equator and Equinox of the date.
!              ir = 3 : Mean Ecliptic of the date and True Equinox.
!              ir = 4 : True Equator and Equinox.
!
!     icd      Index of Coordinates (integer).
!              icd = 1 : Geometric Coordinates.
!              icd = 2 : Apparent Coordinates.
!
!     ipr      Index of Precision (integer).
!              ipr = 1 : High Precision.
!              ipr = 2 : Middle Precision
!              ipr = 3 : Low Precision.
!
! --- Output -----------------------------------------------------------
!
!     c1       Longitude (degree) or Right Ascension (hour) (real*8).
!     c2       Latitude or Declination in degree (real*8).
!     c3       Earth-Moon Distance in kilometer (real*8).
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
! --- Computation  -----------------------------------------------------
!
      select case (ir)
      case (1)
         call MECC (t,c1,c2,c3,icd)
      case (2)
         call MEQC (t,c1,c2,c3,icd)
      case (3)
         call TECC (t,c1,c2,c3,icd,ipr)
      case (4)
         call TEQC (t,c1,c2,c3,icd,ipr)
      endselect
!
      return
      end
!
!
!
      subroutine MECC (t,along,alat,r,icd)
! ----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Object -----------------------------------------------------------
!
!     Mean Ecliptic Coordinates of the Moon.
!     Longitude and Latitude of the Moon, Earth-Moon Distance.
!     Frame : Mean Ecliptic and Equinox of the date.
!
! --- Input ------------------------------------------------------------
!
!     t        Time in Julian century since J2000 (real*8).
!
!     icd      Index of Coordinates (integer).
!              icd = 1 : Geometric Coordinates.
!              icd = 2 : Apparent Coordinates.
!
! --- Output -----------------------------------------------------------
!
!     along    Longitude of the Moon in degree (real*8).
!     alat     Latitude of the Moon in degree (real*8).
!     r        Earth-Moon Distance in kilometer (real*8).
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
      dimension c(3)
!
      common/cv/deg,hr
      common/long/al(5)
      common/lu/pbc(218,3),apbc(218,5,3),per(423,3),aper(423,2,3),
     .          kc(3),kp(3,3)
!
! --- Computation ------------------------------------------------------
!
      do i=1,3
!
         c(i)=0
!
         do j=1,kc(i)
            phi=apbc(j,5,i)
            do j1=4,1,-1
               phi =phi*t+apbc(j,j1,i)
            enddo
            s=sin(phi)
            c(i)=c(i)+s*pbc(j,i)
         enddo
!
         isup=kp(i,1)+kp(i,2)+kp(i,3)
         ct=0
!
         do it=3,1,-1
            inf=isup-kp(i,it)+1
            cc=0
            do j=inf,isup
               phi=aper(j,1,i)+t*aper(j,2,i)
               s=sin(phi)
               cc=cc+s*per(j,i)
            enddo
            ct=cc+t*ct
            isup=inf-1
         enddo
!
         c(i)=c(i)+ct
!
      enddo
!
      phi=al(5)
!
      do i=4,1,-1
         phi=phi*t+al(i)
      enddo
!
      along=c(1)+phi
      alat=c(2)
      r=385000.57d0+c(3)
!
      if (icd.eq.2) then
         phi=(225.d0+477198.9d0*t)*deg
         along=along-0.19524d-3-0.1059d-4*sin(phi)
         r=r+0.708d-1*cos(phi)
         phi=(183.3d0+483202.d0*t)*deg
         alat=alat-0.1754d-4*sin(phi)
      endif
!
      along=mod(along,360.d0)
      if(along.lt.0) along=along+360
!
      return
      end
!
!
!
      subroutine MEQC (t,ra,de,r,icd)
! ----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Object -----------------------------------------------------------
!
!     Mean Equatorial Coordinates of the Moon.
!     Right Ascension and Declination of the Moon, Earth-Moon Distance.
!     Frame : Mean Equator and Equinox of the date.
!
! --- Input ------------------------------------------------------------
!
!     t        Time in Julian century since J2000 (real*8).
!
!     icd      Index of Coordinates (integer).
!              icd = 1 : Geometric Coordinates.
!              icd = 2 : Apparent Coordinates.
!
! --- Output -----------------------------------------------------------
!
!     ra       Right Ascension of the Moon in hour (real*8).
!     de       Declination of the Moon in degree (real*8).
!     r        Earth-Moon Distance in kilometer (real*8).
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
      common/cv/deg,hr
!
! --- Computation ------------------------------------------------------
!
      call OBLI (t,epsi)
!
      epsi=epsi*deg
      ce=cos(epsi)
      se=sin(epsi)
!
      call MECC (t,alg,alt,r,icd)
!
      alg=alg*deg
      alt=alt*deg
      ca=cos(alg)
      sa=sin(alg)
      cb=cos(alt)
      sb=sin(alt)
!
      xq1=cb*ca
      xe2=sa*cb
      xq2=ce*xe2-se*sb
      xq3=se*xe2+ce*sb
!
      ra=atan2(xq2,xq1)
      de=asin(xq3)
      ra=ra/hr
      de=de/deg
      if (ra.lt.0) ra=ra+24
!
      return
      end
!
!
!
      subroutine TECC (t,along,alat,r,icd,ipr)
! ----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Object -----------------------------------------------------------
!
!     True Ecliptic Coordinates of the Moon.
!     Longitude and Latitude of the Moon, Earth-Moon Distance.
!     Frame : Mean Ecliptic of the date and True Equinox.
!
! --- Input ------------------------------------------------------------
!
!     t        Time in Julian century since J2000 (real*8).
!
!     icd      Index of Coordinates (integer).
!              icd = 1 : Geometric Coordinates.
!              icd = 2 : Apparent Coordinates.
!
!     ipr      Index of Precision (integer).
!              ipr = 1 : High Precision.
!              ipr = 2 : Middle Precision.
!              ipr = 3 : Low Precision.
!
! --- Output -----------------------------------------------------------
!
!     along    Longitude of the Moon in degree (real*8).
!     alat     Latitude of the Moon in degree (real*8).
!     r        Earth-Moon Distance in kilometer (real*8).
!
! --- Computation ------------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
      call MECC (t,along,alat,r,icd)
      call NUTA (t,dpsi,deps,ipr)
!
      along=along+dpsi
      along=mod(along,360.d0)
      if(along.lt.0) along=along+360
!
      return
      end
!
!
!
      subroutine TEQC (t,ra,de,r,icd,ipr)
! ----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Object -----------------------------------------------------------
!
!     True Equatorial Coordinates of the Moon.
!     Right Ascension and Declination of the Moon, Earth-Moon Distance.
!     Frame : True Equator and Equinox.
!
! --- Input ------------------------------------------------------------
!
!     t        Time in Julian century since J2000 (real*8).
!
!     icd      Index of Coordinates (integer).
!              icd = 1 : Geometric Coordinates.
!              icd = 2 : Apparent Coordinates.
!
!     ipr      Index of Precision (integer).
!              ipr = 1 : High Precision.
!              ipr = 2 : Middle Precision.
!              ipr = 3 : Low Precision.
!
! --- Output -----------------------------------------------------------
!
!     ra       Right Ascension of the moon in hour (real*8).
!     de       Declination of the Moon in degree (real*8).
!     r        Earth-Moon Distance in kilometer (real*8).
!
! --- Computation-------------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
      common/cv/deg,hr
!
      call OBLI (t,epsi)
      call NUTA (t,dpsi,deps,ipr)
!
      epsi=epsi+deps
      epsi=epsi*deg
      ce=cos(epsi)
      se=sin(epsi)
!
      call MECC (t,alg,alt,r,icd)
!
      alg=alg+dpsi
      alg=alg*deg
      alt=alt*deg
!
      ca=cos(alg)
      sa=sin(alg)
      cb=cos(alt)
      sb=sin(alt)
!
      xq1=cb*ca
      xe2=sa*cb
      xq2=ce*xe2-se*sb
      xq3=se*xe2+ce*sb
!
      ra=atan2(xq2,xq1)
      de=asin(xq3)
      ra=ra/hr
      de=de/deg
      if (ra.lt.0) ra=ra+24
!
      return
      end
!
!
!
      subroutine LUNOSC (t,a,e,gam,pom,gom,alam)
! ----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Object -----------------------------------------------------------
!
!     Osculating Element of the Moon orbit.
!
! --- Input ------------------------------------------------------------
!
!     t        Time in Julian century since J2000 (real*8).
!
! --- Output -----------------------------------------------------------
!
!     a        Semi-major Axis in km (real*8).
!     e        Eccentricity (real*8).
!     gam      Sine of the half inclination (real*8).
!     pom      Longitude of Perigee in degree (real*8).
!     gom      Longitude of Ascending Node in degree (real*8).
!     alam     Mean Longitude of the Moon in degree (real*8).
!
!     Frame : Mean Ecliptic and Equinox of the date.
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
      common/cv/deg,hr
!
      data g/4.0184172d24/
!
! --- Computation ------------------------------------------------------
!
      call DERIV (t,v,u,r,vp,up,rp)
!
      v=v*deg
      u=u*deg
!
      csv=cos(v)
      ssv=sin(v)
      csu=cos(u)
      ssu=sin(u)
!
      x=vp*csu*csu
      y=up*ssv-vp*ssu*csu*csv
      z=up*csv+vp*ssu*csu*ssv
!
      w=sqrt(x*x+y*y+z*z)
      wp=w*deg
      den=sqrt(2*w*(w+x))
      gc=z/den
      gs=y/den
      gam2=gc*gc+gs*gs
      gam=sqrt(gam2)
      gom=atan2(gs,gc)
      gom=gom/deg
      if (gom.lt.0) gom=gom+360.d0
!
      rw=r*wp
      rw2=rw*rw
      a=1.d0/(2.d0/r-(rp*rp+rw2)/g)
      c=ssu/sqrt(1-gam2)
      xp=csu*csv-c*gs
      yp=csu*ssv+c*gc
      zp=rw2*r/g-1
      ws=rp*r*rw/g
      ec=xp*zp+yp*ws
      es=yp*zp-xp*ws
      e2=ec*ec+es*es
      e=sqrt(e2)
      pom=atan2(es,ec)
      pom=pom/deg
      if (pom.lt.0) pom=pom+360.d0
!
      raci=sqrt(1-e2)
      c=ws/(raci*(1+raci))
      rsa=r/a
      cepi=ec+rsa*(xp-es*c)
      sepi=es+rsa*(yp+ec*c)
      epi=atan2(sepi,cepi)
      alam=epi-(ec*sepi-es*cepi)
      alam=alam/deg
      alam=mod(alam,360.d0)
      if (alam.lt.0) alam=alam+360.d0
!
      return
      end
!
!
!
      subroutine OBLI (t,epsi)
!-----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Object -----------------------------------------------------------
!
!     Mean Obliquity.
!
! --- Input ------------------------------------------------------------
!
!     t        Time in Julian century since J2000 (real*8).
!
! --- Output -----------------------------------------------------------
!
!     epsi     Mean Obliquity in degrees (real*8).
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
      dimension eps(5)
!
      data eps/23.43928d0,-0.13d-1,0.d0,0.555d-6,-0.14d-9/
!
! --- Computation ------------------------------------------------------
!
      epsi=eps(5)
!
      do k=4,1,-1
         epsi=epsi*t+eps(k)
      enddo
!
      return
      end
!
!
!
      subroutine NUTA (t,dpsi,deps,ipr)
! ----------------------------------------------------------------------
!
!     Ref : GF-0512
!
! --- Object -----------------------------------------------------------
!
!     Nutation in Longitude and in Obliquity.
!
! --- Input ------------------------------------------------------------
!
!     t        Time in Julian century since J2000 (real*8).
!
!     ipr      Index of Precision (integer).
!              ipr = 1 : High Precision.
!              ipr = 2 : Middle Precision.
!              ipr = 3 : Low Precision.
!
! --- Output -----------------------------------------------------------
!
!     dpsi     Nutation in Longitude in degree (real*8).
!     deps     Nutation in Obliquity in degree (real*8).
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
      dimension dps(13),dep(9),ar(13,3),dps1(5),dep1(2),kn(3)
!
      common/cv/deg,hr
!
! --- Initialization ---------------------------------------------------
!
      data kn/13,2,1/
!
      data dps/
     . -4.77767d-3,-0.36631d-3,-0.6317d-4, 0.5728d-4,
     .   0.3961d-4,  0.1978d-4,-0.1436d-4,-0.1072d-4, -0.836d-5,
     .    0.603d-5,  -0.439d-5,  0.358d-5, 0.342d-5/
!
      data dps1/
     . -0.484d-5,3*0.d0,-0.9d-7/
!
      data dep/
     . 2.55625d-3,0.15933d-3,0.2714d-4,-0.2486d-4,
     . 2*0.d0,0.622d-5,0.556d-5,0.358d-5/
!
      data dep1/
     . 0.25d-6,-0.9d-7/
!
      data ar/
     .    125.0446d0,   200.9329d0,    76.6333d0,    250.0891d0,
     .    357.5291d0,   134.9634d0,   198.4620d0,    311.5888d0,
     .    211.5967d0,
     .    203.4038d0,   259.2630d0,     75.8883d0,   301.6699d0,
     .  -1934.1362d0, 72001.5397d0, 962535.7627d0, -3868.2724d0,
     .  35999.0503d0,477198.8676d0, 108000.5899d0,964469.8989d0,
     .1439734.6303d0, 36002.4894d0,-413335.3554d0, 73935.6758d0,
     . 485336.8951d0,
     .  21.d-4,6.d-4,-27.d-4, 42.d-4,  -2.d-4,90.d-4,5.d-4,-47.d-4,
     .  63.d-4,8.d-4,123.d-4,-15.d-4,-117.d-4/
!
      dpsi=0
      deps=0
!
! --- Computation ------------------------------------------------------
!
      do i=1,kn(ipr)
!
         phi=ar(i,3)
!
         do j=2,1,-1
            phi=phi*t+ar(i,j)
         enddo
!
         phi=phi*deg
         c=dps(i)
         if (i.eq.1.or.i.eq.5) c=c+t*dps1(i)
         dpsi=dpsi+c*sin(phi)
!
         if (i.ne.5.and.i.ne.6.and.i.le.9) then
            c=dep(i)
            if (i.eq.1.or.i.eq.2) c=c+t*dep1(i)
            deps=deps+c*cos(phi)
         endif
!
      enddo
!
      return
      end
!
!
!
      subroutine DERIV (t,v,u,r,vp,up,rp)
!-----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Object -----------------------------------------------------------
!
!     Mean Coordinates of the Moon and Derivatives.
!
! --- Input ------------------------------------------------------------
!
!     t        Time in Julian century since J2000 (real*8).
!
! --- Output -----------------------------------------------------------
!
!     v        Longitude in degree (real*8).
!     u        Latitude in degree (real*8).
!     r        Distance in km (real*8).
!
!     vp       Longitude Derivative in degree/cy (real*8).
!     up       Latitude Derivative in degree/cy (real*8).
!     rp       Distance Derivative in km/cy (real*8).
!
!     Frame : Mean Ecliptic and Equinox of the date.
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
      dimension c(3),dc(3),dal(4)
!
      common/lu/pbc(218,3),apbc(218,5,3),per(423,3),aper(423,2,3),
     .          kc(3),kp(3,3)
      common/long/al(5)
!
      data dal/481267.881d0,-0.26536d-2,0.5568d-5,-0.6136d-7/
!
! --- Computation ------------------------------------------------------
!
      do i=1,3
!
         c(i)=0
         dc(i)=0
!
         do j=1,kc(i)
            phi=apbc(j,5,i)
            do j1=4,1,-1
               phi=phi*t+apbc(j,j1,i)
            enddo
            dphi=4*apbc(j,5,i)
            do j1=4,2,-1
               dphi=dphi*t+(j1-1)*apbc(j,j1,i)
            enddo
            c(i)=c(i)+pbc(j,i)*sin(phi)
            dc(i)=dc(i)+pbc(j,i)*dphi*cos(phi)
         enddo
!
         isup=kp(i,1)+kp(i,2)+kp(i,3)
         dct=0
         ct=0
!
         do it=3,1,-1
            inf=isup-kp(i,it)+1
            cc=0
            dcc=0
            do j=inf,isup
               phi=aper(j,1,i)+t*aper(j,2,i)
               cc=cc+per(j,i)*sin(phi)
               dcc=dcc+per(j,i)*aper(j,2,i)*cos(phi)
            enddo
            ct=cc+t*ct
            dct=dcc+t*dct
            isup=inf-1
         enddo
!
         c(i)=c(i)+ct
         dc(i)=dc(i)+dct
!
      enddo
!
      phi=al(5)
      do i=4,1,-1
         phi=phi*t+al(i)
      enddo
!
      v=c(1)+phi
      u=c(2)
      r=385000.57d0+c(3)
      v=mod(v,360.d0)
      if (v.lt.0) v=v+360.d0
!
      phi=dal(4)
!
      do i=3,1,-1
         phi=phi*t+dal(i)
      enddo
!
      vp=dc(1)+phi
      up=dc(2)
      rp=dc(3)
!
      return
      end
!
!
!
      subroutine TITRE (iprint,nul,ipr,icd,ir,idt,idf)
!-----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Object -----------------------------------------------------------
!
!     Printing the titles of the results.
!
! --- Input ------------------------------------------------------------
!
!     iprint   Index of Printing (integer).
!              iprint = 1 : Screen.
!              iprint = 2 : File.
!
!     nul      Logical unit number of the results file (integer).
!
!     ipr      Index of Precision (integer).
!              ipr = 1 : High Precision.
!              ipr = 2 : Middle Precision
!              ipr = 3 : Low Precision.
!
!     icd      Index of Coordinates (integer).
!              icd = 1 : Geometric Coordinates.
!              icd = 2 : Apparent Coordinates.
!              icd = 3 : Osculating Elements.
!
!     ir       Index of Frame (integer).
!              ir = 1 : Mean Ecliptic Frame.
!              ir = 2 : Mean Equator Frame.
!              ir = 3 : True Ecliptic Frame.
!              ir = 4 : True Equator Frame.
!
!     idt      Index of Time (integer).
!              idt = 1 : Dynamical Time.
!              idt = 2 : Universal Time.
!
!     idf      Index of Dates (integer).
!              idt = 1 : Julian Dates.
!              idt = 2 : Calendar Dates.
!
! --- ------------------------------------------------------------------
!
      select case (iprint)
!
      case (1)
!
         call EFFECR
         write (*,1000)
!
         select case (icd)
         case (1)
            write (*,1001) 
!V98        write (*,1001) ipr
         case (2)
            write (*,1002) 
!V98        write (*,1002) ipr
         case (3)
            write (*,1003) 
!V98        write (*,1003) ipr
            write (*,2001)
            if (idf.eq.1) write (*,5001)
            if (idf.eq.2) write (*,5002)
            return
         endselect
!
         select case (ir)
         case (1)
            write (*,2001)
         case (2)
            write (*,2002)
         case (3)
            write (*,2003)
         case (4)
            write (*,2004)
         endselect
!
         ncas=(idf-1)*4+(idt-1)*2+2-mod(ir,2)
         select case (ncas)
         case (1)
            write (*,3001)
         case (2)
            write (*,3002)
         case (3)
            write (*,3003)
         case (4)
            write (*,3004)
         case (5)
            write (*,4001)
         case (6)
            write (*,4002)
         case (7)
            write (*,4003)
         case (8)
            write (*,4004)
         endselect
!
      case (2)
!
         select case (icd)
         case (1)
            write (nul,1001) 
!V98        write (nul,1001) ipr
         case (2)
            write (nul,1002) 
!V98        write (nul,1002) ipr
         case (3)
            write (nul,1003) 
!V98        write (nul,1003) ipr
            write (nul,2001)
            if (idf.eq.1) write (nul,5001)
            if (idf.eq.2) write (nul,5002)
            return
         endselect
!
         select case (ir)
         case (1)
            write (nul,2001)
         case (2)
            write (nul,2002)
         case (3)
            write (nul,2003)
         case (4)
            write (nul,2004)
         endselect
!
         ncas=(idf-1)*4+(idt-1)*2+2-mod(ir,2)
         select case (ncas)
         case (1)
            write (nul,3001)
         case (2)
            write (nul,3002)
         case (3)
            write (nul,3003)
         case (4)
            write (nul,3004)
         case (5)
            write (nul,4001)
         case (6)
            write (nul,4002)
         case (7)
            write (nul,4003)
         case (8)
            write (nul,4004)
         endselect
!
      endselect
!
      return
!
! --- Formats ----------------------------------------------------------
!
1000  format (2x,'JANUS GF9812',15x,'MOON from -4000 to 8000',16x,
     .        'M&J CHAPRONT'/2x,78('-'))
!
1001  format (/2x,t31,'GEOMETRIC COORDINATES')
1002  format (/2x,t32,'APPARENT COORDINATES ')
1003  format (/2x,t32,'OSCULATING ELEMENTS')
!V98  1001  format (/2x,t23,'GEOMETRIC COORDINATES (precision',i2,')')
!V98  1002  format (/2x,t24,'APPARENT COORDINATES (precision',i2,')')
!V98  1003  format (/2x,t20,'OSCULATING ELEMENTS (precision',i2,')')
!
2001  format (2x,t22,'Mean Ecliptic and Equinox of the date'/)
2002  format (2x,t23,'Mean Equator and Equinox of the date'/)
2003  format (2x,t20,'Mean Ecliptic of the date and True Equinox'/)
2004  format (2x,t27,'True Equator and True Equinox'/)
!
3001  format (2x,t12,'Julian Date (ET)',
     .        8x,'Longitude',5x,'Latitude',7x,'Distance')
3002  format (2x,t12,'Julian Date (ET)',
     .        7x,'Right Asc.',4x,'Declination',5x,'Distance')
3003  format (2x,t6,'Julian Date (UT)',
     .        8x,'Longitude',5x,'Latitude',7x,'Distance',
     .        5x,'ET-UT')
3004  format (2x,t6,'Julian Date (UT)',
     .        7x,'Right Asc.',4x,'Declination',5x,'Distance',
     .        5x,'ET-UT')
!
4001  format (2x,t10,1x,'D',2x,'M',3x,'Yr',5x,'H (ET)',
     .        5x,'Longitude',5x,'Latitude',7x,'Distance')
4002  format (2x,t10,1x,'D',2x,'M',3x,'Yr',5x,'H (ET)',
     .        4x,'Right Asc.',4x,'Declination',5x,'Distance')
4003  format (2x,t4,1x,'D',2x,'M',3x,'Yr',5x,'H (UT)',
     .        5x,'Longitude',5x,'Latitude',7x,'Distance',
     .        5x,'ET-UT')
4004  format (2x,t4,1x,'D',2x,'M',3x,'Yr',5x,'H (UT)',
     .        4x,'Right Asc.',4x,'Declination',5x,'Distance',
     .        5x,'ET-UT')
!
5001  format (4x,'Julian Date (ET)',
     .        10x,'a',8x,'e',5x,'Gamma',
     .        3x,'L_Peri',3x,'L_Node',4x,'Lambda')
5002  format (3x,'D',2x,'M',3x,'Yr',5x,'H (ET)',
     .        7x,'a',8x,'e',5x,'Gamma',
     .        3x,'L_Peri',3x,'L_Node',4x,'Lambda')
!
      end
!
!
!
      subroutine IMPC (d,nul,along,alat,r,ipr,dt,idt,idf)
!-----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Object -----------------------------------------------------------
!
!     Printing the Ecliptic Coordinates.
!
! --- Input ------------------------------------------------------------
!
!     d        Julian Date (real*8).
!     nul      Logical unit number of the results file (integer).
!     ipr      Index of Precision (integer).
!     dt       Difference ET-UT in seconds (real*8).
!     idt      Index of Time (integer).
!     idf      Index of Date (integer).
!
! --- Output -----------------------------------------------------------
!
!     along    Longitude in degree (real*8).
!     alat     Latitude in degree (real*8).
!     r        Earth-Moon Distance in km (real*8).
!
! --- ------------------------------------------------------------------
!
      implicit double precision (a-h,o-z)
      character*22 adat
      character*50 acd
!
      id=int(dt)
!
      select case (idf)
      case (1)
         write (adat,1001) d
      case (2)
         call DJDATE (d,3,ij,im,ia,nh,nm,ns)
         write (adat,1002) ij,im,ia,nh,nm,ns
      endselect
!
      select case (idt)
      case (1)
         select case (ipr)
         case (1)
            write (acd,2001) along,alat,r
            if (acd(5:6).eq.' .') acd(5:6)='0.'
            if (acd(18:19).eq.' .') acd(18:19)='0.'
            if (acd(18:19).eq.'-.') acd(17:19)='-0.'
         case (2)
            rp=r/10.d0
            write (acd,2002) along,alat,rp
            if (acd(6:7).eq.' .') acd(6:7)='0.'
            if (acd(19:20).eq.' .') acd(19:20)='0.'
            if (acd(19:20).eq.'-.') acd(18:20)='-0.'
            acd(37:37)='0'
         case (3)
            rp=r/100.d0
            write (acd,2003) along,alat,rp
            if (acd(7:8).eq.' .') acd(7:8)='0.'
            if (acd(20:21).eq.' .') acd(20:21)='0.'
            if (acd(20:21).eq.'-.') acd(19:21)='-0.'
            acd(36:37)='00'
         endselect
         write (*,2004) adat,acd
         write (nul,2004) adat,acd
      case (2)
         select case (ipr)
         case (1)
            write (acd,3001) along,alat,r,id
            if (acd(5:6).eq.' .') acd(5:6)='0.'
            if (acd(18:19).eq.' .') acd(18:19)='0.'
            if (acd(18:19).eq.'-.') acd(17:19)='-0.'
         case (2)
            rp=r/10.d0
            write (acd,3002) along,alat,rp,id
            if (acd(6:7).eq.' .') acd(6:7)='0.'
            if (acd(19:20).eq.' .') acd(19:20)='0.'
            if (acd(19:20).eq.'-.') acd(18:20)='-0.'
            acd(37:37)='0'
         case (3)
            rp=r/100.d0
            write (acd,3003) along,alat,rp,id
            if (acd(7:8).eq.' .') acd(7:8)='0.'
            if (acd(20:21).eq.' .') acd(20:21)='0.'
            if (acd(20:21).eq.'-.') acd(19:21)='-0.'
            acd(36:37)='00'
         endselect
         write (*,3004) adat,acd
         write (nul,3004) adat,acd
      endselect
!
      return
!
! --- Formats ----------------------------------------------------------
!
1001  format (4x,f13.5,5x)
1002  format (i2.2,1x,i2.2,1x,i5,2x,i2.2,'h',i2.2,'m',i2.2,'s')
!
2001  format (f10.4,'d',2x,f10.4,'d',2x,f11.1,' km')
2002  format (f10.3,'d',2x,f10.3,'d',2x,f11.0,' km')
2003  format (f10.2,'d',2x,f10.2,'d',2x,f10.0,'  km')
2004  format (t10,a22,2x,a40)
!
3001  format (f10.4,'d',2x,f10.4,'d',2x,f11.1,' km',2x,i6,' s')
3002  format (f10.3,'d',2x,f10.3,'d',2x,f11.0,' km',2x,i6,' s')
3003  format (f10.2,'d',2x,f10.2,'d',2x,f10.0,'  km',2x,i6,' s')
3004  format (t4,a22,2x,a50)
!
      end
!
!
!
      subroutine IMPQ (d,nul,ra,de,r,ipr,dt,idt,idf)
!-----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Object -----------------------------------------------------------
!
!     Printing the Ecliptorial Coordinates.
!
! --- Input ------------------------------------------------------------
!
!     d        Julian Date (real*8).
!     nul      Logical unit number of the results file (integer).
!     ipr      Index of Precision (integer).
!     dt       Difference ET-UT in seconds (real*8).
!     idt      Index of Time (integer).
!     idf      Index of Date (integer).
!
! --- Output -----------------------------------------------------------
!
!     ra       Right Ascension in hour (real*8).
!     de       Declination un degree (real*8).
!     r        Earth-Moon Distance in km (real*8).
!
! --- ------------------------------------------------------------------
!
      implicit double precision (a-h,o-z)
      character*22 adat
      character*50 acd
!
      id=int(dt)
!
      select case (idf)
      case (1)
         write (adat,1001) d
      case (2)
         call DJDATE (d,3,ij,im,ia,nh,nm,ns)
         write (adat,1002) ij,im,ia,nh,nm,ns
      endselect
!
      select case (idt)
      case (1)
         select case (ipr)
         case (1)
            write (acd,2001) ra,de,r
            if (acd(4:5).eq.' .') acd(4:5)='0.'
            if (acd(18:19).eq.' .') acd(18:19)='0.'
            if (acd(18:19).eq.'-.') acd(17:19)='-0.'
         case (2)
            rp=r/10.d0
            write (acd,2002) ra,de,rp
            if (acd(5:6).eq.' .') acd(5:6)='0.'
            if (acd(19:20).eq.' .') acd(19:20)='0.'
            if (acd(19:20).eq.'-.') acd(18:20)='-0.'
            acd(37:37)='0'
         case (3)
            rp=r/100.d0
            write (acd,2003) ra,de,rp
            if (acd(6:7).eq.' .') acd(6:7)='0.'
            if (acd(20:21).eq.' .') acd(20:21)='0.'
            if (acd(20:21).eq.'-.') acd(19:21)='-0.'
            acd(36:37)='00'
         endselect
         write (*,2004) adat,acd
         write (nul,2004) adat,acd
      case (2)
         select case (ipr)
         case (1)
            write (acd,3001) ra,de,r,id
            if (acd(4:5).eq.' .') acd(4:5)='0.'
            if (acd(18:19).eq.' .') acd(18:19)='0.'
            if (acd(18:19).eq.'-.') acd(17:19)='-0.'
         case (2)
            rp=r/10.d0
            write (acd,3002) ra,de,rp,id
            if (acd(5:6).eq.' .') acd(5:6)='0.'
            if (acd(19:20).eq.' .') acd(19:20)='0.'
            if (acd(19:20).eq.'-.') acd(18:20)='-0.'
            acd(37:37)='0'
         case (3)
            rp=r/100.d0
            write (acd,3003) ra,de,rp,id
            if (acd(5:6).eq.' .') acd(5:6)='0.'
            if (acd(19:20).eq.' .') acd(19:20)='0.'
            if (acd(19:20).eq.'-.') acd(18:20)='-0.'
            acd(36:37)='00'
         endselect
         write (*,3004) adat,acd
         write (nul,3004) adat,acd
      endselect
!
      return
!
! --- Formats ----------------------------------------------------------
!
1001  format (4x,f13.5,5x)
1002  format (i2.2,1x,i2.2,1x,i5,2x,i2.2,'h',i2.2,'m',i2.2,'s')
!
2001  format (f10.5,'h',2x,f10.4,'d',2x,f11.1,' km')
2002  format (f10.4,'h',2x,f10.3,'d',2x,f11.0,' km')
2003  format (f10.3,'h',2x,f10.2,'d',2x,f10.0,'  km')
2004  format (t10,a22,2x,a40)
!
3001  format (f10.5,'h',2x,f10.4,'d',2x,f11.1,' km',2x,i6,' s')
3002  format (f10.4,'h',2x,f10.3,'d',2x,f11.0,' km',2x,i6,' s')
3003  format (f10.3,'h',2x,f10.2,'d',2x,f10.0,'  km',2x,i6,' s')
3004  format (t4,a22,2x,a50)
!
        end
!
!
!
      subroutine IMPS (d,nul,a,e,gam,pom,gom,alam,ipr,idf)
!-----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Object -----------------------------------------------------------
!
!     Printing the Osculating Elements.
!
! --- Input ------------------------------------------------------------
!
!     d        Julian Date (real*8).
!     nul      Logical unit number of the results file (integer).
!     ipr      Index of Precision (integer).
!     idf      Index of Date (integer).
!
! --- Output -----------------------------------------------------------
!
!     a        Semi-major Axis in km (real*8).
!     e        Eccentricity (real*8).
!     gam      Sine of the half inclination (real*8).
!     pom      Longitude of Perigee in degree (real*8).
!     gom      Longitude of Ascending Node in degree (real*8).
!     alam     Mean Longitude of the Moon in degree (real*8).
!
! --- ------------------------------------------------------------------
!
      implicit double precision (a-h,o-z)
      character*22 adat
      character*54 acd
!
      select case (idf)
      case (1)
         write (adat,1001) d
      case (2)
         call DJDATE (d,3,ij,im,ia,nh,nm,ns)
         write (adat,1002) ij,im,ia,nh,nm,ns
      endselect
!
      select case (ipr)
      case (1)
         write (acd,2001) a,e,gam,pom,gom,alam
         acd(7:7)='k'
         acd(10:10)='0'
         acd(18:18)='0'
         if (acd(30:31).eq.' .') acd(30:31)='0.'
         if (acd(38:39).eq.' .') acd(38:39)='0.'
         if (acd(47:48).eq.' .') acd(47:48)='0.'
      case (2)
         ap=a/100.d0
         write (acd,2002) ap,e,gam,pom,gom,alam
         acd(5:7)='00k'
         acd(11:11)='0'
         acd(19:19)='0'
         if (acd(31:32).eq.' .') acd(31:32)='0.'
         if (acd(39:40).eq.' .') acd(39:40)='0.'
         if (acd(49:50).eq.' .') acd(49:50)='0.'
      case (3)
         ap=a/100.d0
         write (acd,2003) ap,e,gam,pom,gom,alam
         acd(5:7)='00k'
         acd(12:12)='0'
         acd(20:20)='0'
         acd(34:34)='d'
         if (acd(40:41).eq.' .') acd(40:41)='0.'
         if (acd(50:51).eq.' .') acd(50:51)='0.'
      endselect
!
      write (*,3000) adat,acd
      write (nul,3000) adat,acd
!
      return
!
! --- Formats ----------------------------------------------------------
!
1001  format (4x,f13.5,5x)
1002  format (i2.2,1x,i2.2,1x,i5,2x,i2.2,'h',i2.2,'m',i2.2,'s')
!
2001  format (f7.0,'m',1x,f7.5,1x,f8.6,2x,
     .        f6.2,'d',1x,f7.3,'d',1x,f8.4,'d')
2002  format (f5.0,'  m',1x,f7.4,1x,f8.5,2x,
     .        f6.1,'d',1x,f7.2,'d',1x,f8.2,'d')
2003  format (f5.0,'  m',1x,f7.3,1x,f8.4,2x,
     .        f7.0,1x,f7.1,'d',1x,f8.2,'d')
!
3000  format (2x,a22,2x,a54)
!
        end
!
!
!
      subroutine DCDJ (cd,d,ierr,datcal)
!-----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Object -----------------------------------------------------------
!
!     Change a Calendar Date in a Julian Date.
!
! --- Input ------------------------------------------------------------
!
!     cd       Calendar date +/-yyyymmdd.hhmmss (real*8).
!
! --- Output -----------------------------------------------------------
!
!     d        Julian Date (real*8).
!
!     ierr     Index of Error (integer).
!              ierr = 0 : No Error
!              ierr = 1 : Date between 05 Oct 1582 and 15 Oct 1582.
!              ierr = 2 : Not a valid date.
!
!     datcal   Calendar Date if ierr=0 or 1 (character*30).
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
      character*4 cmo(12)
      character*16 car
      character*30 datcal
!
! --- Initialization ---------------------------------------------------
!
      data cmo/
     .'Janv','Febr','Mars','Apr.','May','June',
     .'July','Aug.','Sept','Oct.','Nov.','Dec.'/
!
      d=0.d0
      ierr=0
      datcal=' '
!
! --- Decode the Calendar Date -----------------------------------------
!
      write (car,1000) cd
      if (car(10:10).ne.'.') then
         ierr=2
         return
      endif
!
      read (car(1:5),1001) iy
      read (car(6:7),1002) im
      read (car(8:9),1002) id
      read (car(11:12),1002) ih
      read (car(13:14),1002) imn
      read (car(15:16),1002) is
      read (car(1:9),1003) iw
!
! --- Control of the date parameters -----------------------------------
!
      if (im.lt.1.or.im.gt.12.or.
     .    id.gt.31.or.ih.gt.24.or.imn.gt.59.or.is.gt.59) then
          ierr=2
          return
      endif
!
      write (datcal,2000) id,cmo(im),iy,ih,imn,is
!
! --- Compute the Julian date ------------------------------------------
!
      if (im.le.2) then
         ims=im+13
         iyp=iy-1
      else
         ims=im+1
         iyp=iy
      endif
!
      if (iw.ge.15821015) then
         ia=iyp/100
         ib=2-ia+ia/4
      else
         if (iw.ge.15821005) then
            ierr=1
            return
         endif
         ib=0
      endif
!
      w=365.25d0*iyp
      if (iyp.lt.0) w=w-0.75d0
!
      d=int(w)+int(30.60001d0*ims)+id+ib+ih/24.d0+imn/1440.d0
     . +is/86400.d0+1720994.5d0
!
      return
!
! --- Formats ----------------------------------------------------------
!
1000  format(f16.6)
1001  format(i5)
1002  format(i2)
1003  format(i9)
2000  format(i2.2,1x,a4,1x,i5,3x,i2.2,'h',i2.2,'m',i2.2,'s')
!
        end
!
!
!
      subroutine DJDATE (tdj,icode,jour,mois,kan,nh,nm,ns)
!-----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Object -----------------------------------------------------------
!
!     Change a Julian Date in a Calendar Date.
!
! --- Input ------------------------------------------------------------
!
!     tdj         Julian Date (real*8).
!
!     icode       Index of conversion code (integer)
!                 icode=0 : Conversion tdj -> day-month-year & 00-00-00.
!                 icode=1 : Conversion tdj -> day-month-year & nh-00-00.
!                 icode=2 : Conversion tdj -> day-month-year & nh-nm-00.
!                 icode=3 : Conversion tdj -> day-month-year & nh-nm-ns.
!
! --- Output -----------------------------------------------------------
!
!     jour        Day number (integer).
!     mois        Month number (integer).
!     kan         Year number (integer).
!
!     nh          Hour number (integer).
!     nm          Minutes number (integer).
!     ns          Seconds (number integer).
!
! --- Declarations -----------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
! --- Initialization ---------------------------------------------------
!
      jour=0
      mois=0
      kan=0
      nh=0
      nm=0
      ns=0
!
      if (tdj.lt.0.d0.or.icode.lt.0.or.icode.gt.3) return
      if (icode.eq.0) dt=0.0d0
      if (icode.eq.1) dt=0.5d0/24.d0
      if (icode.eq.2) dt=0.5d0/1440.d0
      if (icode.eq.3) dt=0.5d0/86400.d0
!
! --- Conversion -------------------------------------------------------
!
      t=tdj+dt+0.5d0
      z=int(t)
      f=t-z
!
      if (z.lt.2299161.d0) then
         a=z
      else
         alpha=int((z-1867216.25d0)/36524.25d0)
         a=z+1.d0+alpha-int(alpha/4.d0)
      endif
!
      b=a+1524.d0
      c=int((b-122.1d0)/365.25d0)
      d=int(365.25d0*c)
      e=int((b-d)/30.6001d0)
!
      jour=b-d-int(30.6001d0*e)
!
      if (e.lt.13.5d0) then
         mois=e-1.d0
      else
         mois=e-13.d0
      endif
!
      if (mois.lt.2.5d0) then
         kan=c-4715.d0
      else
         kan=c-4716.d0
      endif
!
      f=f*24.d0
      nh=f
      f=(f-nh)*60.d0
      nm=f
      f=(f-nm)*60.d0
      ns=f
!
      return
      end
!
!
!
      subroutine UT (dj,dt,ierr,limut,tlimut)
! ----------------------------------------------------------------------
!
!     Ref : GF-0212
!
! --- Object -----------------------------------------------------------
!
!     Difference ET-UT.
!
! --- Input ------------------------------------------------------------
!
!     dj       Julian Date (real*8).
!
! --- Output -----------------------------------------------------------
!
!     dt       ET-UT in seconds (real*8).
!
!     ierr     Index of Error (integer).
!              ierr = 0  : No Error.
!              ierr = 1  : Date greater than the UT limit.
!
!     limut    Year of the UT limit (integer).
!     tlimut   Julian date of the UT limit (integer).
!
! --- Remarks ----------------------------------------------------------
!
!     The Julian Date 'dj' may be expressed indefferently in ET or UT.
!
!     For the dates before 1980, the values are derived from Stevenson 
!     and Morrison (1984) : before 1600, the values of ET-UT are given
!     by formulae, after 1600, they are stored in the table 'idtsm' 
!     (every 5 years before 1780 and every year after 1780).
!
!     For the dates since 1980 onwards, we use yearly values derived 
!     from observed TAI-UT1 (IERS); they are also stored in the table 
!     'idtsm'. 
!
!     After the year given by the parameter 'limit', there is no   
!     evaluation of ET-UT (dt=0). 
!
! --- Yearly Update ----------------------------------------------------
!
!     - Add 1 to the size of the table 'idtsm' (parameter idim),
!     - Add the new value of ET-TU (in second) in the table 'idtsm'.
!
! --- ------------------------------------------------------------------
!
      implicit double precision (a-h,o-z)
!
! *** ******************************************************************
!     Last release : 2005, ET-TU=64 seconds (idim=262)
!
      parameter (idim=262)
      dimension idtsm(idim)
      data idtsm/
     .  80,78,75,73,70,68,65,52,43,36,30,26,20,16,10,
     .  06,01,-2,-3,-4,-3,-3,-2,-1,00,01,01,03,03,05,
     .  05,06,08,09,09,11,11,
     .  11,11,11,11,12,12,12,12,12,12,12,11,11,11,11,
     .  10,10,10,09,09,09,09,08,08,08,08,08,08,08,08,
     .  09,09,09,09,09,09,09,09,09,09,08,08,08,07,07,
     .  06,06,06,05,05,04,04,03,03,03,03,03,03,03,03,
     .  03,03,04,04,04,04,04,05,05,05,05,05,05,06,06,
     .  06,06,06,06,06,06,06,05,05,04,04,03,01,00,00,
     .  -1,-2,-3,-4,-4,-5,-6,-6,-6,-6,-6,-6,-6,-6,-7,
     .  -7,-7,-7,-6,-7,-7,-7,-7,-7,-7,-7,-6,-5,-4,-3,
     .  -2,-1,01,02,03,05,06,07,09,10,11,13,14,16,17,
     .  18,19,20,21,21,22,22,23,23,23,24,24,24,24,24,
     .  24,24,24,24,24,24,24,24,24,24,25,25,26,26,27,
     .  27,28,28,29,29,30,30,30,31,31,31,32,32,33,33,
     .  34,34,34,35,36,37,37,38,39,40,41,42,43,44,45,
     .  46,48,49,50,51,51,52,53,54,54,55,55,56,56,57,
     .  58,58,59,60,61,62,63,63,64,64,64,64,64,64,64/
!
! *** ******************************************************************
!
      dimension f1(3),f2(3)
      data f1/2177.d0,495.d0,42.4d0/
      data f2/102.d0,100.d0,23.6d0/
!
      data d2000/2451545.d0/
      data d1780/2371191.5d0/
      data sjul/36525.d0/
!
      data t948/-10.519658d0/
      data t1600/-3.999932d0/
      data t1780/-2.199959d0/
!
      ierr=0
      dt=0.d0
!
      ia=idim-36
      limut=1780+ia
      tlimut=d1780+ia*sjul/100.d0
!
      if (dj.ge.tlimut) then
         ierr=1
         return
      endif
!
      t=(dj-d2000)/sjul
!
      if ((t-t948).lt.0) then
         dt=(f1(3)*t+f1(2))*t+f1(1)
      else
         if ((t-t1600).lt.0) then
            dt=(f2(3)*t+f2(2))*t+f2(1)
         else
            if ((t-t1780).lt.0) then
               i=int((t-t1600)*20)+1
            else
               i=int((t-t1780)*100)+37
            endif
            dt=idtsm(i)
         endif
      endif
!
      return
      end
!
!
!
      subroutine EFFECR
! ----------------------------------------------------------------------
!
!     Ref : GF-9812
!
! --- Object -----------------------------------------------------------
!
!     Clear the screen and position the cursor.
!     Depending of the operating system : 'cls' or 'clear'.
!
! ----------------------------------------------------------------------
!
      call SYSTEM ('cls')
!!!   call SYSTEM ('clear')
      write (*,'(/)')
!
      return
      end
      